!
! MMM - see description by author P.B.Stetson
!       Fortran 90 version by FH

!------------------------------------------------------------------------------

Module MMModule

! global parameters
      Integer, parameter :: MAXIT = 30
      Integer, parameter :: MINARRAY = 5

contains

    Subroutine MMM(narray, array, HiBad, mean, median, mode, sigma, skew)

! Input:
!     array ... sorted array of the input values 
!     HiBad ... upper cut-off
!
! Output:
!     narray ... No. of used values
!     mean   ... mean
!     median ... median
!     mode   ... mode
!     sigma  ... standard deviation
!     skew   ... skewness

      Integer :: narray
      Real :: array(narray), HiBad, mean, median, mode, sigma, skew

!#######################################################################
!
!      SUBROUTINE  MMM (SKY, NSKY, HIBAD, SKYMN, SKYMED,
!     .     SKYMOD, SIGMA, SKEW)
!
!=======================================================================
!
!               Official DAO version:  1988 July 1
!
! This version of MMM (modified by PBS 1984.IV.10ff) assumes that
! the sky brightnesses in the one-dimensional array SKY are already
! sorted on entering this routine, and that pixels outside the "bad"
! limits have already been eliminated.
!
! This particular version of MMM also takes cognizance of the fact that,
! pixels falling below the LOBAD threshold already having been
! eliminated, the contaminated sky pixels values overwhelmingly display
! POSITIVE departures from the true value.
!
! If for some reason it is impossible to obtain the mode of the sky
! distribution, this will be flagged by setting SIGMA = -1.0.
!
! Arguments
!
!     SKY (INPUT) is a real vector containing actual sorted sky values.
!    NSKY (INPUT) is the number of defined elements in SKY.
!  SKYMOD (OUTPUT) is the estimated mode of the sky values.
!   SIGMA (OUTPUT) is the computed standard deviation of the peak in
!         the sky histogram.
!    SKEW (OUTPUT) is the computed skewness of the peak in the sky
!         histogram.
!
!=======================================================================
!

!      Integer, parameter :: MAXIT = 30
!      Integer, parameter :: MINARRAY = 20

      Double Precision ::  sum, sumsq
      REAL :: cut, cut1, cut2, delta, mid
      REAL :: r, x, center, side
      Integer :: i, minimm, maximm, niter, istep, jstep
      Logical :: redo
!
! Section 1
!
!      narray = size(array)
      if( narray < MINARRAY ) goto 9900
!
      mid=0.5*(array((narray + 1)/2) + array(narray/2 + 1))
!
! SKYMID is the median value for the whole ensemble of sky pixels.
! Notice that if NSKY is odd, then (NSKY+1)/2 and (NSKY/2)+1 will be the
! same number, whereas if NSKY is even, they will be different numbers.
! This same trick will be used again later.
!
! Initialize the variables for accumulating the mean and standard
! deviation, and initialize the rejection limits.
!
      sum = 0.0
      sumsq = 0.0
      cut1 = min(mid - array(1), array(narray) - mid, HiBad - mid)
!
! For the first pass we will consider only pixels in a symmetric
! interval of brightness values about the median value.  This exploits
! the assumption that all the bad pixels are already rejected from the
! lower end of the brightness range.
!
      cut2 = mid + cut1
      cut1 = mid - cut1; !write(*,*) cut1,cut2
!
      minimm = 0
      maximm = narray
      Do i  = 1,narray
         if( array(i) < cut1 )then
            minimm = i
         else if( array(i) > cut2)then
            exit
         else
           delta = array(i) - mid
           sum = sum + delta
           sumsq = sumsq + delta**2
           maximm = i
         endif
      enddo; !write(*,*) minimm, maximm
!
! Henceforth in this subroutine, MINIMM will point to the highest value
! rejected at the lower end of the vector, and MAXIMM will point to the
! highest value accepted at the upper end of the vector.
! MAXIMM-MINIMM is the number of pixels within the acceptance range.
!
! Compute mean and sigma (from the first pass).
!
      median = 0.5*(array((minimm + maximm + 1)/2) &
                 +  array((minimm + maximm)/2+1))
      mean = sum / (maximm - minimm)
      sigma = sqrt(sumsq/(maximm - minimm) - mean**2)
      mean = mean + mid
!
! The middle sky value, SKYMID, was subtracted off up above and added
! back in down here to reduce the truncation error in the computation
! of SIGMA.
! Note that this definition of SIGMA is incorrect by a factor of
! SQRT [NSKY/(NSKY-1.)], but for all but pathological cases (where none
! of this can be trusted anyway), it's close enough.
!
      mode = mean
      if(median < mean) mode = 3.0*median - 2.0*mean
!
! If the mean is less than the mode, that means the contamination is
! slight, and the mean value is what we really want.  Note that this
! introduces a slight bias toward underestimating the sky when
! the scatter in the sky is caused by random fluctuations rather than
! by contamination, but I think this bias is negligible compared to the
! problem of contamination.
!
! Section 2
!
! Rejection and recomputation loop:
!
      niter = 0
      Do 
        niter = niter + 1
        if( niter > MAXIT .or. (maximm - minimm) < MINARRAY ) goto 9900
!
! Compute Chauvenet rejection criterion.
!
        r = log10( float(maximm - minimm))
        r = max(2.0, (-0.1042*r + 1.1695)*r + 0.8895)
!
! Compute rejection limits (symmetric about the current mode).
!
        cut = r*sigma + 0.5*abs(mean - mode)
        cut = max(1.5 , cut)
        cut1 = mode - cut
        cut2 = mode + cut
!
! Recompute mean and sigma by adding and/or subtracting sky values
! at both ends of the interval of acceptable values.
!
! At each end of the interval, ISTEP will show the direction we have to
! step through the vector to go from the old partition to the new one.
! Pixels are added or subtracted depending upon whether the limit is
! moving toward or away from the mode.
!
        redo = .FALSE.
!
! Is CUT1 above or below the minimum currently-accepted value?
!
        istep = int(sign(1.0001,cut1 - array(minimm + 1)))
        jstep = (istep + 1)/2
!
! If ISTEP = +1, JSTEP = 1.  If ISTEP = -1, JSTEP=0.  If ISTEP = +1,
! then we know that at least one pixel must be deleted at the low end.
!
        if( istep > 0)  goto 2120
 2100   if( istep < 0 .and. minimm <= 0 ) goto 2150
!
! Quit when SKY(MINIMM) < CUT1 <= SKY(MINIMM+1)
!
          if( array(minimm) <= cut1 .and. array(minimm + 1) >= cut1) goto 2150
!
! If ISTEP is positive, subtract out the sky value at MINIMM+1; if
! ISTEP is negative, add in the sky value at MINIMM.
!
2120      continue
          delta = array(minimm + jstep) - mid
          sum = sum - istep*delta
          sumsq = sumsq - istep*delta**2
          minimm = minimm + istep
          redo = .TRUE.                                ! A change has occured
        goto 2100
2150    continue
!write(*,*) istep
!
! Is CUT2 above or below the current maximum?
!
        istep = int(sign(1.0001, cut2 - array(maximm)))
        jstep = (istep + 1)/2
!
! If ISTEP = +1, JSTEP = 1.  If ISTEP = -1, JSTEP=0.  If ISTEP = -1,
! then we know that we must subtract at least one pixel from the high
! end.
!
        if(istep < 0) goto 2220
 2200   if( istep > 0 .AND. maximm >= narray ) goto 2250
!
! Quit when SKY(MAXIMM) <= CUT2 < SKY(MAXIMM+1)
!
          if( array(maximm) <= cut2 .and. array(maximm + 1) >= cut2) goto 2250
!
! If ISTEP is positive, add in the sky value at MAXIMM+1; if ISTEP is
! negative, subtract off the sky value at MAXIMM.
!
 2220     continue
          delta = array(maximm + jstep) - mid
          sum = sum + istep*delta
          sumsq = sumsq + istep*delta**2
          maximm = maximm + istep
          redo = .TRUE.                                 ! A change has occured
        goto 2200
!
 2250   continue
!write(*,*) '2>',istep
!
! Compute mean and sigma (from this pass).
!
        mean = sum/(maximm - minimm)
        sigma = sqrt(sumsq/(maximm - minimm) - mean**2)
        mean = mean + mid
!
! Obtain the median.  To first approximation, the median would be the
! value of the sky in the middle pixel in the sorted data (if the
! total number is odd) or the mean of the two pixels straddling
! the middle (if the total number of pixels is even).
!
!     SKYMED=0.5*(SKY((MINIMM+MAXIMM+1)/2)+SKY((MINIMM+MAXIMM)/2+1))
!
! However, this is not good enough.  If you look at the estimator for
! the mode, you will note that a tiny change in the list of sky pixels,
! just sufficient to alter the median value of the sky brightness by
! one unit, will change the estimator of the mode by three units.  We
! really want something more robust than this.  As a first attempt
! at a more robust median estimator, I propose to estimate the median
! of the distribution by the mean of the central five percent of sky
! values.  This involves considerable care to make sure you get
! a perfectly symmetric sample of pixels about the median, whether
! there is an even or an odd number of pixels within the acceptance
! interval.
!
        median = 0.0
        x = 0.0
        center = real(minimm + 1 + maximm)/2.0
        side = real(nint(0.05*(maximm - minimm)))/2. + 0.25
!
        Do i = nint(center - side), nint(center + side)
          median = median + array(i)
          x = x + 1.0
        enddo
!
        median = median/x
        mode = mean
        if(median < mean) mode = 3.0*median - 2.0*mean
!
! If the mean is less than the mode, that means the contamination is
! slight, and the mean value is what we really want.  Note that this
! introduces a slight bias toward underestimating the sky when
! the scatter in the sky is caused by random fluctuations rather than
! by contamination, but I think this bias is negligible compared to the
! problem of contamination.
!
! If the limits have not yet stopped moving, try again.
!
!        write(*,*) redo
        if( .not. redo ) exit
      enddo
!
!-----------------------------------------------------------------------
!
! Normal return.
!
      skew = (mean - mode)/max(1.0, sigma)
      narray = maximm - minimm
!      write(*,*) narray
      Return
!
!-----------------------------------------------------------------------
!
! An error condition has been detected.
!
 9900 sigma = -1.0
      skew = 0.0
!      write(*,*) sigma
      Return
!
      End Subroutine
!
!----------------------------------------------------------------------------

End 
